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Abstract 

i> 

Q\ \ This paper discusses the possibility of constructing time-independent solutions 

to the collisionless Boltzmann equation which depend on quantities other than 
global isolating integrals such as energy and angular momentum. The key point 
is that, at least in principle, a self-consistent equilibrium can be constructed 
from any set of time-independent phase space building blocks which, when 
^r) • combined, generate the mass distribution associated with an assumed time- 

independent potential. This approach provides a way to justify Schwarzschild's 
(1979) method for the numerical construction of self-consistent equilibria with 
■ arbitrary time- independent potentials, generalising thereby an approach devel- 

oped by Vandervoort (1984) for integrable potentials. As a simple illustration, 
Schwarzschild's method is reformulated to allow for a straightforward compu- 
Q\ ■ tation of equilibria which depend only on one or two global integrals and no 

other quantities, as is reasonable, e.g., for modeling axisymmetric configura- 
tions characterised by a nonintegrable potential. 
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O 1 1 MOTIVATION 

CO 

Conventional wisdom holds that galaxies in or near equilibrium can be modeled as time- 
independent solutions to the collisionless Boltzmann equation. In this context, the modeling 
of galaxies would seem to break logically into two reasonably distinct components, namely 
(i) constructing time-independent solutions to the collisionless Boltzmann equation and 
then (ii) determining whether said solutions are linearly stable and otherwise viable as 
reasonable models of what one actually sees. The principal focus of this paper is primarily 
on the former component, the construction of time-independent solutions, although the 
concluding section will comment on issues related to viability. 

Given an assumed equilibrium mass distribution po, and an associated potential $o 
generated by the gravitational Poisson equation, 

V 2 $ = 4vrGp , (1) 

there is a globally conserved quantity (isolating integral) E reflecting time translational 
invariance. If the configuration is time-independent in an inertial frame, this quantity is 
the ordinary energy E = ^v 2 + <I>o (here, and henceforth, units have been chosen so that the 
stellar mass m = 1). If instead the configuration is time- independent in a suitably chosen 
rotating frame, E is the Jacobi integral. If <J>o manifests other continuous symmetries as 
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well (e.g., spherical or axial symmetry), Noether's Theorem (cf. Arnold 1989) implies that 
there will be one or more additional globally conserved quantities, say {ij}. If there exist 
three independent global integrals, motion in the potential $0 is integrable. If fewer than 
three independent global integrals exist, the motion is nonintegrable. 

Jeans' Theorem implies that any function fo(E, {h}) that reproduces the assumed po, 
i.e., for which 

Po = J tfVo(£,Ui}), (2) 

yields a self-consistent equilibrium (cf. Binney and Tremaine 1987). Indeed, many work- 
ers have gone further and assumed that all time-independent equilibria must depend on 
such global isolating integrals. If, however, one demands that every /q be a function only 
of some set of global isolating integrals, there is an obvious critique which can be leveled 
towards numerical model-building based on Schwarzschild's (1979) method. This method, 
which involves selecting ensembles of orbit segments that self-consistently reproduce the 
mass density associated with an imposed potential, says nothing a priori about any global 
integrals; and, as such, for generic potentials nothing intrinsic to the method imposes ex- 
plicitly the demand that the equilibrium /q be a function of one or more global isolating 
integrals. 

For the special case of integrable potentials, e.g., spherical configurations or triaxial equi- 
libria characterised by Staeckel (1890) potentials, there is a direct, albeit not completely 
trivial, connection between Schwarzschild's method and global integrals. As discussed by 
Vandervoort (1984) in a slightly different language, if orbits are constrained by three in- 
dependent global integrals, say Ji(r, v), l2(r,v), and fixing the values of these 
integrals as l2,o, and I^q determines completely a collection of one or more multiply 
periodic orbits in the three-dimensional configuration space, each of which must in principle 
be included in an equilibrium model with the proper relative weight.^] 

More precisely, specifying a triplet {Ji,o, ^2,0 > ^3,0} defines a phase space density 

5/i, ,/2,o,/3,o( r ' v ) « 5 D [h(r,\) - Ji j0 ] 5 D [h(r,\) - I 2 ,o] S D [I 3 (r,\) - J 3j0 ] (3) 
and a corresponding configuration space density 

™/i,o,/ 2 ,o,/3,o( r ) = / d 3 vg h0jl20>l30 (r,v) 

=111 dhdhdh |lxfy ff/i ' o ' /2 ' o ' /3 ' o(r ' v) ' (4) 

which, when evolved into the future using the Hamilton equations associated with <3?o, re- 
mains unchanged. The assumption that the desired po is generated from some fo(h, h,^) 

2 Consider, for example, a spherical system with global integrals E, J 2 , and J z , and focus on 
motion in the equatorial plane, i.e., J x — J y = 0. Here there are an infinite number of possible 
orbits, characterised by initial conditions corresponding to all possible points (r, ip) located in an 
annulus with inner and outer radii fixed by E and J 2 . The weighting implicit in eq. (3) means 
that all values of ip should be treated equally, but that the relative weighting of different r's must 
reflect the fact that, because of conservation of angular momentum, orbits spend different amounts 
of time at different radii. As a practical matter, however, this subtlety is arguably unimportant, 
and it may suffice computationally to consider a single orbit: If the radial and azimuthal periods 
are incommensurate, any (x, y) yields an orbit that densely fills the annulus with the proper weight; 
and, even if the periods are not incommensurate, unless they are in a relatively low order resonance 
the orbit will generally fill a region which, to the level of accuracy associated with one's configuration 
space discretisation, is essentially dense in the annulus. 
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depending only on global integrals means that the equilibrium distribution must be con- 
structed as a superposition of solutions of the form given by eq. (3). In this context, the 
proper construction of an equilibrium model using Schwarzschild's method thus entails three 
stages, namely (i) selecting all the orbits entering into each gi 1}0 ,i 2 ,o,h,o( r i v ) with the weights 
implicit in eq. (3), (ii) performing the integration of eq. (4) to extract nj 1 j 2 ,j 3 (r), which 
incorporates the fact that each point in configuration space is weighted in proportion to the 
amount of time orbits spend there, and then (hi) choosing a superposition of ni 1 j 2 j 3 (r)'s 
which yields the imposed po- 

At least in the specific setting described by Vandervoort (1984), this interpretation 
breaks down for the case of nonintegrable potentials, where one cannot identify three global 
integrals, or, more generally, whenever one relaxes the demand that /o be realisable as a 
function of global isolating integrals. However, as described more carefully later on in this 
paper, one can still capture the essential aspect of Vandervoort's analysis, namely that 
appropriately identified orbit segments yield the natural time-independent building blocks 
in terms of which to construct a self-consistent model. 

Nonintegrable potentials generically admit two different classes of orbits, namely regu- 
lar and chaotic. Regular orbits in a nonintegrable potential behave qualitatively like orbits 
in an integrable potential in that they are multiply periodic and, even more importantly, 
are restricted to a three- (or lower-) dimensional hypersurface in the six-dimensional phase 
space. It follows that, even though they do not admit three global integrals, they must be 
constrained (cf. Lichtenberg and Lieberman 1992) by what are sometimes termed "local in- 
tegrals."^] For this reason, regular orbits in nonintegrable potentials define time-independent 
building blocks in the same sense as do orbits in integrable potentials. Chaotic orbits are 
very different in that they are intrinsically aperiodic and densely fill phase space regions 
that are necessarily higher than three-dimensional. However, chaotic phase space regions 
can still be characterised by invariant distributions which, when evolved into the future, 
remain unchanged; and one can use these invariant distributions as an additional set of 
building blocks when constructing self-consistent equilibria. 

This alternative "orbital" interpretation of the building blocks for self-consistent equilib- 
ria is important for at least two reasons. (1) It is possible to generate analytically exact two- 
integral models for axisymmetric configurations which are characterised by nonintegrable 
potentials (cf. Hunter and Qian 1993), including potentials where motion in the meridional 
plane is chaotic. However, reproducing these models numerically using Schwarzschild's 
method must entail sampling a collection of time-independent building blocks more com- 
plex than those associated with integrable potentials. (2) Because global integrals are 
associated with continuous symmetries, one might expect generically that, for genuinely 
three-dimensional potentials, there will not exist any global integral aside from the energy 
(or the Jacobi integral). If, however, one demands that the equilibrium distribution be a 
function only of the energy E, so that fo = fo(E), one concludes (cf. Binney and Tremaine 
1987) that the mass distribution po must be spherical. (This is the analogue of the well 
known theorem from stellar structure that all static perfect fluid stars are spherical.) In 

3 A good example of a local integral is the so-called third integral associated in some cases with a 
nonintegrable, time-independent, axisymmetric potential, where the only global integrals are energy 
E and rotational angular momentum J z . How regular orbits in a nonintegrable potential differ from 
orbits in an integrable potential is clearly stated on p. 49 of Lichtenberg and Lieberman (1992): 
"Since the regular trajectories depend discontinuously on initial conditions, their presence does not 
imply the existence of an isolating integral (global invariant) or symmetry of the system. However, 
regular trajectories, when they exist, represent exact invariants of the motion." 
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other words, every nonrotating triaxial equilibrium must depend on something other than 
the energy, either additional global integrals, as for Staeckel potentials, or "local" integrals, 
as is implicit in the construction of cuspy triaxial models by Merritt and Fridman (1996) 
or Siopis (1997). 

Section 2 discusses more carefully the basic building blocks of a self-consistent equilib- 
rium, allowing explicitly for the possibility of nonintegrable potentials that admit chaotic 
orbits. Section 3 then illustrates how, for the simple case of nonintegrable equilibria de- 
pending only on one or two global isolating integrals (and nothing else), Schwarzschild's 
method can be reformulated in terms of an appropriate set of time-independent building 
blocks. Section 4 concludes by describing straightforward generalisations to construct equi- 
libria that do not depend simply on global integrals and then discussing the issue of whether 
such "more complex" equilibria are physically viable. 

2 INVARIANT DISTRIBUTIONS AND SELF-CONSISTENT EQUILIBRIA 

It is often asserted glibly that "any equilibrium solution fo to the collisionless Boltzmann 
equation must be given as a function of the time-independent integrals of the motion asso- 
ciated with the potential $o generated self-consistently from fo" In point of fact, however, 
this statement is an oversimplification and requires some careful thought. Should one de- 
mand, as is often assumed, at least tacitly, that fo depends only on the global isolating 
integrals, such as energy E (or the Jacobi integral for a rotating configuration) and angu- 
lar momentum J z , or can one instead allow for equilibria fo that depend on the "local" 
isolating integrals which, in a generic nonintegrable potential, make regular orbits regular, 
i.e., restrict them to a lower-dimensional phase space hypersurface (cf. Lichtenberg and 
Lieberman 1992)? Could one, for example, try to construct models which assign regular 
and chaotic orbits on the same E-J z hypersurface different weights, or must one sample 
each constant E-J z hypersurface uniformly? 

Arguably, the only crucial point is that a time-independent solution to the collisionless 
Boltzmann equation must be constructed from time-independent building blocks, so as to 
ensure that, if initial data be evolved into the future along the characteristics associated 
with the self-consistent potential, the form of fo will remain unchanged. The easiest way 
to do this, both conceptually and practically, is to demand that fo be given as a function 
of one or more global isolating integrals, say E and /. The obvious point is that any 
fo(E, I) which implies the proper mass density po, and hence the proper potential wm 
yield a time-independent solution since, by assumption, both E and I are time-independent 
constants of the motion. In other words, dE/dt = dl/dt = implies dfo/dt = 0. This is 
of course the standard way of showing that time-independent conserved quantities can be 
used to construct a self-consistent equilibrium. 

However, there is another, more "microscopic," viewpoint. Specifically, viewed in terms 
of the orbits associated with the equilibrium (i.e., characteristics associated with the Boltz- 
mann equation), this construction works because such an fo implies that the phase space 
number density is constant on hypersurfaces of constant E and /. This constancy means 
that the orbit ensemble that generates fo must involve a uniform, i.e., microcanonical, sam- 
pling of each constant E — I phase space hypersurface, but Hamilton's equations for motion 
in a fixed $o imply that such a population is invariant under time translation. 

It would seem clear from this latter viewpoint that choosing fo to be a function only of 
global isolating integrals is not necessary, at least in principle. A priori, a self-consistent 
equilibrium fo can be constructed from any collection of time-independent building blocks 
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which successfully reproduces the assumed potential $o- The key point, then, as stressed, 
e.g., by Ott (1993), is that, assuming the validity of the Ergodic Theorem, for flows in a fixed 
time-independent potential any orbit is ergodic in an appropriately interpreted subspace, 
so that a microcanonical population of the appropriate subspace yields a time-independent 
building block. 

Regular orbits are multiply periodic and, as such, are characterised (in an orbit-averaged 
sense) by a density that is time- independent, so that they can be treated individually as 
time-independent building blocks, with a density p proportional to the time that the orbit 
spends in the neighborhood of each point. In this sense, regular orbits in a nonintegrable 
potential can be exploited in the same way as the orbits in an integrable potential, even 
though it is seemingly impossible to identify explicitly the forms of the "local integrals," 
and even though regularity is not attributable directly to a continuous symmetry. 

Chaotic orbits are not periodic, so that this naive argument does not hold. However, it 
would still seem possible to identify an appropriate set of time-independent chaotic building 
blocks. For a fixed value of E (and any other global integral), the constant E (or E-I) phase 
space hypersurface divides naturally into regular and chaotic regions. The chaotic region 
divides in turn into one or more subregions which are connected in the sense that an orbit 
starting in any part of a subregion will eventually pass arbitrarily close to every other part of 
that subregion. The important point then is that a uniform, i.e., microcanonical, population 
of any connected chaotic domain defines a time-independent building block. Why this 
should be so is easy to understand: Time translation using Hamilton's equations moves 
each phase space point in the chaotic domain to another point in the same domain, but the 
only initial distribution invariant under time translation using Hamilton's equations is a 
constant density distribution. Integrating this phase space building block over the allowed 
range of velocities yields a configuration space density which, as for the regular orbits, is 
proportional to the amount of time that a representative chaotic orbit in the domain spends 
in the neighborhood of each point r. 

Assuming the validity of the Ergodic Theorem for individual connected chaotic domains, 
it is relatively simple to generate such invariant distributions numerically. Specifically, one 
knows that, when evolved into the future, a generic ensemble of initial conditions located 
anywhere in the domain will evidence a coarse-grained approach towards the invariant 
microcanonical distribution. In fact, this has been confirmed by numerical experiments (cf. 
Kandrup and Mahon 1994, Mahon, Abernathy, Bradley, and Kandrup 1995, Merritt and 
Valluri 1996) which have shown that, for chaotic flows in a variety of different potentials, 
reduced distribution functions (like /(r) or /(v)) exhibit an apparent exponential in time 
approach towards an invariant reduced distribution on a time scale r that correlates with 
the value of the largest Lyapunov exponent. 

Viewed in this fashion, regular and chaotic orbits can be used to define time-independent 
building blocks in exactly the same way, the only difference being that chaotic building 
blocks are intrinsically higher-dimensional. When evolved into the future, initial conditions 
corresponding to regular and chaotic orbits both yield trajectories which, in an asymptotic, 
late time limit, converge towards time- independent invariant distributions. 

The crucial point in all of this is that, in principle, a library comprised of all possible 
invariant distributions, both those corresponding to individual regular orbits and those 
corresponding to individual connected chaotic phase space regions, should constitute a 
complete set of building blocks in terms of which to construct self-consistent models of 
a galaxy with the specified potential $o- In the real world, one cannot construct such 
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a library, which would contain an infinite number of building blocks. However, one can 
construct a large, but finite, library and then sample that library in an attempt to select 
appropriate combinations that reproduce a suitably discretised version of the assumed <!>o- 
This is the essence of what Schwarzschild's (1979) method can, and should, do when applied 
to a generic nonintegrable potential that admits both regular and chaotic orbits. It is also 
evident that, in principle, nothing stops one from trying to construct equilibria that contain 
only regular (or perhaps only chaotic) orbits, although one might imagine that it would be 
very hard to reproduce a smooth potential $0 with a collection of orbits that systematically 
avoids significant phase space regions or, especially, probes the phase space in an exceedingly 
irregular fashion. 

For the case of a generic rotating, axisymmetric equilibrium, there are two global iso- 
lating integrals, namely E and J z , associated respectively via Noether's Theorem with 
symmetries with respect to time translations and rotations about the z-axis. The potential 
may in fact be integrable, so that there are three global isolating integrals, but this is not 
necessarily the case. In general, a rotating, axisymmetric, time-independent potential will 
be nonintegrable and admit both chaotic and regular orbits, each of which defines time- 
independent building blocks. It follows that, if one allows for local integrals, one can, at 
least in principle, try to construct equilibria that do not sample constant E-J z surfaces uni- 
formly. One could, e.g., try to construct models which exclude all chaotic orbits. Analytic 
approaches to constructing self-consistent axisymmetric equilibria, as developed, e.g., by 
Hunter and Qian (1994), neglect this possibility altogether and focus exclusively on solu- 
tions for which /o = fo(E, J z ), so that the phase space density is constant on hypersurfaces 
of constant E and J z . Whether this is well motivated physically, or whether this is only a 
useful analytic simplification, is not completely clear at the present time. 

It should be stressed that great care must be taken in identifying the invariant distribu- 
tions associated with (ensembles of) chaotic orbits. As has long been known from numerical 
investigations of simple maps (cf. Lieberman and Lichtenberg 1972), the presence of cantori 
(cf. Aubry and Andre 1978, Mather 1982) or an Arnold (1964) web allows for the possibility 
of chaotic near-invariant distributions which, albeit not strictly time- independent, can, at 
least in the absence of any perturbations, behave for very long times as if they were essen- 
tially time-independent distributions. As discussed, e.g., in Mahon, Abernathy, Bradley, 
and Kandrup (1995), the crucial point here that cantori and Arnold webs serve as partial 
obstructions which, although they cannot completely block motion between different phase 
space regions, can significantly impede phase space transport. It follows that, even if a 
single chaotic region is connected, it may appear partitioned into disjoint regions even over 
relatively long time scales. 

This phenomenon is problematic. A putative self-consistent equilibrium generated with 
a near-, rather than true, invariant distribution cannot be a true self-consistent equilibrium. 
On sufficiently long time scales, the orbital population associated with the distribution will 
change, occasioning changes in the mass distribution, the potential, and so forth. One 
might nevertheless want to argue that, if the time scale associated with this phenomenon 
is sufficiently long, this very slow effect will be irrelevant astronomically, so that one can 
speak of nearly self-consistent equilibria that can exist for times much longer than the age 
of the Universe, tu- However, this argument is probably wrong. Real astronomical systems 
involve iV-body realisations of self-consistent equilibria which, heuristically, are presumed 
to behave like smooth three-dimensional Hamiltonian systems perturbed by friction and 
noise. However, numerical experiments involving perturbations of motion in a fixed poten- 
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tial indicate (Habib, Kandrup, and Mahon 1996, 1997) that even very weak friction and 
noise can dramatically accelerate phase space transport through cantori or along an Arnold 
web, occasioning significant changes in an initial near- invariant distribution on compara- 
tively short time scales. Trying to estimate the longevity of a near- invariant distribution 
without allowing for the effects of perturbations is unquestionably a very bad idea. 

3 A VARIANT OF SCHWARZSCHILD'S METHOD FOR ONE- AND TWO- 
INTEGRAL DISTRIBUTIONS 

The objective of this Section is to reformulate Schwarzschild's method in terms of the natural 
set of building blocks so as to permit the construction of equilibrium models /o which are 
assumed to depend on one or two global integrals and to exhibit no additional dependence 
on any nonclassical local integrals. This is a straightforward generalisation of an approach 
proposed by Vandervoort (1984) for the construction of three-integral equilibria. 

Start by specifying a time-independent potential <3?o( r )) an d hence a configuration space 
density /9o(r),0 which admits (say) two constants of the motion, E and I, where E is the 
particle energy (or, perhaps, the Jacobi integral) and I is some other isolating integral, the 
form of which is assumed to be known explicitly.^ By assumption, the desired equilibrium 
(or equilibria) Jq must be given exactly as a function /q = fo(E,I). The object, therefore, 
is to construct a discretised approximation to a smooth /o of this form which reproduces 
the assumed p(r) self-consistently. This can be done in two stages, viz: 

1. First grid E-I space into a collection of cells and, for each pair {Ei-Ij}, write down 
the invariant distribution gij(r, v) on the constant Ei-Ij hypersurface. Use these <?«'s 
to derive reduced configuration space densities ny(r). 

2. Then construct the desired numerical approximation to /o as a sum of contributions 
from the different invariant distributions g^j , with the relative weights of the different 
gij's fixed by the requirement of self-consistency for the configuration space density. 

This construction proceeds without explicit reference to individual orbits and, as such, 
provides no insight into the orbital building blocks entering into the equilibrium. If this 
be perceived as a serious lacuna, the natural tack numerically is to consider separately the 
different constant E-I hypersurfaces and, on each hypersurface, to construct ensembles of 
orbit segments that reproduce self-consistently the invariant <7«'s. 

3.1 Construction of the invariant distribution for fixed E and I 

A uniform population of the phase space hypersurface of fixed E{ and Ij corresponds to 
an invariant distribution of the form 

g(E u Ij) = gy(r,v) = K S D [Ei - E(r,v)] S D [Ij - 7(r, v)], (5) 

where 5d denotes a Dirac delta, and the quantities E and I are viewed explicitly as functions 
of the phase space coordinates. The quantity K, is a constant, whose value is fixed by the 

4 Although it is the density, rather than the potential, that astronomers are wont to specify, 
it is more natural conceptually to view 4>o as the fundamental object, since it is the Hamiltonian 
associated with 4>o that defines the time-invariant building blocks. 

5 Reformulating the following for equilibria admitting only one isolating integral is completely 
straightforward. If the equilibrium admits three independent integrals, it is integrable and can be 
addressed using the approach described by Vandervoort (1984). 
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normalisation 

J d 3 x J d 3 vg tJ (r,v) = 1, (6) 

where the integral extends over the allowed phase space regions. In other words, the in- 
variant distribution corresponds to a normalised microcanonical population of the constant 
Ei-Ij hypersurface. 

Given such a gij(r, v), it is straightforward to integrate over the velocity dependence to 
extract a reduced configuration space density nij(r). Because E and I are known functions 
of r and v, one can choose to view any two of the phase space coordinates, say v y and v z , 
as functions of E, I, and the remaining four phase space coordinates.^] However, the Dirac 
deltas in eq. (5) make dv y and dv z integrations trivial, so that one can immediately write 
down analytically a reduced 

9ij(x,y,z,v x ) = K J dv y dv z 5 D [Ei - E(r,v)] 5 D [Ij - I(r,v)]. (7) 

It follows that the configuration space density, 

mj(x,y,z)= dv x gij(x,y,z,v x ), (8) 



associated with each constant {Ei,Ij} pair is given as a simple quadrature. 

In general, it may be impossible to perform the integral in eq. (8) analytically. This, 
however, is not a serious difficulty. Even if known analytically, the n^s must eventually be 
approximated by a set of values on a configuration space grid so as to facilitate a comparison 
between the imposed density po(r) and the inferred density n(r) constructed from the 
invariant n^s. 

3.2 Construction of fo(E,I) from the invariant distributions 

In the continuum limit, one knows that the true equilibrium distribution 

/o(r,v)=| J dEdIA(E,I) g E ,i(r,v), (9) 

where qej is the invariant distribution for fixed E and I, viewed as a function of r and 
v, and A(E, I) is an expansion coefficient, which gives the relative weights of the different 
values of E and / entering into /o- The discretised construction thus involves 

/ (r,v)=££A^(r,v). (10) 

* 3 

The proper choice of weights Aij derives from the demand of self-consistency: Given /o, 
one can define a density 

n(r) = J d 3 vf (r,v) (11) 

which, when discretised, becomes 

n{x,y,z) = ^2^2A ij n ij (x,y,z) (12) 

* 3 



6 The choice of Cartesian coordinates, implicit in the following, is only for specificity: as far as 
this algorithm is concerned, the coordinate system is completely irrelevant. 
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in terms of the unknown expansion coefficients A^. However, demanding that this n(x, y, z) 
correspond as closely as possible to the density 

po(x,y,z) = ^ V 2 $ (13) 

associated with the assumed potential <&o then enables one to determine the "best" values for 
the Aij's. This construction is very much analogous to the ordinary Schwarzschild method, 
save only that the basic building blocks are now the reduced invariant distributions fiij, 
rather than individual orbits. 

3.3 Orbital building blocks for the invariant distributions 

One way in which to obtain insights into the orbital building blocks of a self-consistent 
model constructed using this algorithm is by proceeding numerically to construct ensem- 
bles of orbit segments which reproduce self-consistently the invariant distributions gij. In 
general, will contain contributions from both regular and chaotic orbits, each of which 
is characterised separately by its own invariant distribution. The easiest way to construct 
gij is probably to (1) obtain an invariant distribution for the chaotic orbits and then (2) 
augment this by another (sub) distribution comprised of segments of regular orbits, the 
latter so chosen that the combination of regular and chaotic orbits yields a satisfactory 
approximation to the true invariant distribution. 

The invariant distribution is approximated numerically by binning the six-dimen-sional 
phase space into a collection of six-dimensional hypercubes, and then assigning occupation 
numbers to the different hypercubes which are proportional to the time that orbits sampling 
the true invariant distribution reside in each cell. (This is justified by the Ergodic Theorem 
[cf. Lichtenberg and Lieberman 1992].) In the continuum limit, the invariant distribu- 
tion corresponds to a uniform population on a four-dimensional phase space hypersurface. 
Given a discretisation of the phase space coordinates, the invariant distribution corresponds 
instead to a four-dimensional shell in the six-dimensional phase space. 

The invariant (sub-)distribution of chaotic orbits is especially easy to compute if, as 
is often the case, for fixed Ei and Ij the entire chaotic region is connected in the sense 
that, eventually, every chaotic orbit will pass arbitrarily close to every point in the chaotic 
region. (For simplicity, ignore the tiny measure of chaotic orbits trapped permanently inside 
invariant KAM tori.) All that one need do is specify a (more or less arbitrary) ensemble of 
initial conditions, each corresponding to a chaotic orbit, evolve each initial condition into 
the future, and wait until the evolved ensemble approaches an invariant distribution, i.e., 
a uniform sampling of the chaotic portions of the Ei-Ij hypersurface (cf. Kandrup and 
Mahon 1994, Mahon, Abernathy, Bradley, and Kandrup 1995). 

To expedite the calculation, it is useful to evolve the initial conditions in the presence 
of very weak amplitude friction and noise, sufficiently weak that the values of E and / arc 
almost constant (cf. Habib, Kandrup, and Mahon 1996, 1997). The advantage of intro- 
ducing weak friction and noise is that such small perturbations can dramatically accelerate 
the overall approach towards a true invariant distribution by facilitating extrinsic diffusion 
through cantori and/or along an Arnold web (cf. Lichtenberg and Lieberman 1992). If one 
does not either (a) integrate for a very long time and/or (b) allow for such perturbing influ- 
ences, one faces the problem that the initial ensemble may evolve towards a near-invariant 
distribution which, albeit not strictly time-independent, only changes significantly on a very 
long time scale. 

The contribution of different regular orbits to the invariant distribution can be generated 
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using an analogue of the original Schwarzschild method. Specify a large number of regular 
initial conditions and integrate each into the future to generate a library of regular orbits. 
Then use a linear programming algorithm, or some variant thereof, to select a weighted 
ensemble of regular orbits which, when combined with the chaotic (sub)distribution, yields 
a satisfactory approximation to the true invariant distribution on the constant Ei-Ij hyper- 
surface. 

This construction of invariant distributions gij, and the corresponding densities n^-, 
is admittedly tedious numerically (albeit presumably straightforward), since it involves 
repeating Schwarzschild 's method for each pair {Ei,Ij}. However, it is arguably a crucial 
step in obtaining a proper understanding of the orbital structure associated with the self- 
consistent model since, as discussed already, one knows that the g^s are the proper building 
blocks in terms of which to construct an equilibrium fo(E, I). 

3.4 A simple two-dimensional example 

Consider as a pedagogical example the case of two-dimensional gravity, this corresponding 
physically to a collection of infinite rods aligned along the z-axis, and suppose that the 
configuration is rotating uniformly about the z-axis with angular velocity CI. It then follows 
that, in the rotating frame, the configuration is characterised by a potential ^(x,y) and a 
surface density a(x, y) related by 

V 2 ^(x,y) = 4irGa(x,y). (15) 

Suppose then that there is only one global isolating integral, namely the Jacobi integral 
E, which, in terms of phase space coordinates defined in the rotating frame, takes the form 

E = \vl + \v 2 y + V(x, y) - Ul\x 2 + y 2 ) = ±v 2 x + \v 2 y + * eff (x, y). (16) 

To the extent that one demands that any equilibrium /o associated with this mass distri- 
bution be a function only of the isolating integral E, the fundamental building block is the 
microcanonical phase space density on a constant Jacobi integral hypersurface, which, for 
any Ei, takes the form 

g(Ei) = (fc(r,v) = KS D [Ei - E(r,v)]. (17) 
As will be evident from below, the normalisation constant K, can be written as 

C =^) <18> 

in terms of V(Ei), the area of the configuration space region with ^f e ff < Ei- The reduced 
configuration space density rij associated with this gi satisfies 

rii(x,y)=IC J J dv x dv y 5 D [Ei - E(r,v)\, (19) 

where the integrals extend over the values of v x and v y that are allowed energetically, i.e., for 
which Ei > ^ e ff- The dv y integration can be performed trivially by implementing the Dirac 
delta, allowing for nonzero contributions at two values, namely v y = ±^2[E~i — ^ e ff) ~ v x = 
a. It follows that, for those regions in configuration space for which fy e ff(x,y) < Ei, 

n(x,y) = 2/C f -^t= . (20) 
J -a \/a z — V* 
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The remaining integral can then be performed trivially, leading to a reduced configuration 
space density on the constant Ei hypersurface of the form 



ni(x,y) = 2TrK,Q[Ei - * e// (x,y)] 



1 



e[Ei-V eff (x,y)], 



(21) 



V(Ei) 



where Q(z) = 1 for z > and = otherwise. 

It follows from eq. (21) that, independent of the specific form of the potential ^f(x,y), 
the total configuration space surface density, given as a sum of contributions on different 
constant Jacobi integral hypersurfaces, must be of the form 



where the A^'s give the relative weights of the different hypersurfaces. The demand that 
the n(x, y) of eq. (22) agree as closely as possible with the a(x, y) associated with ^ may 
then be used to identify the "best" values of the A^s. 

4 DISCUSSION 

There are a number of different ways in which the algorithm described in the preceding 
Section can be generalised to permit the construction of more complex equilibria, which 
do not depend simply on the global integrals E and /. For fixed values of E and /, it is 
straightforward to locate the general locations of (at least the large) chaotic regions and, by 
evolving arbitrary ensembles of initial conditions located in these regions into the future, 
it is easy to derive a numerical approximation to the invariant distribution associated with 
each of these chaotic regions. Given these invariant distributions, one can then integrate 
over velocities to extract the chaotic contribution ji^-(x, y, z) to total density riij(x,y, z) 
associated with any pair Ei and Ij. Subtracting from the full n-ij then yields the 
regular contribution n^-(x, y, z) to the density. However, given a knowledge of n\j and n^- 
separately, one can then attempt to construct models which assign different relative weights 
to the regular and chaotic portions of the Ei-Ij hypersurface, thus allowing one to test the 
prejudice of some workers that self-consistent models should contain few, if any, chaotic 
orbits. 

Similarly, one can identify those portions of the Ei-Ij hypersurface that correspond to 
different types of regular orbits, e.g., boxes and tubes, and compute their relative densities, 
say n\j and n\j, which can in turn be used as separate building blocks. In particular, given 
such a collection one can try to construct models which associate different relative weights 
to boxy and/or tuby and/or chaotic orbits, and, to the extent that such models can be 
constructed, one can investigate whether the different phase space densities /o have obvious 
observational signatures which could be compared with real astronomical data. Is there, 
e.g., some natural signature which, when observed in real galaxies, can be interpreted as 
evidence that /q contains a significant measure of chaotic orbits? 

In principle, one can continue this process of refinement more or less ad infinitum, 
identifying increasing numbers of time-independent building blocks associated with different 
regular orbits, although one's freedom to deal with chaotic orbits is limited by the fact that 
there is only one natural notion of a time-independent invariant distribution. However, it is 
not clear that such a procedure is well motivated physically. At least heuristically, it would 
seem that building an equilibrium by "picking and choosing" amongst individual orbits in 
a strongly nonintegrable potential with different values of local isolating integrals is akin to 




(22) 



11 



selecting orbits in an integrable potential which yield a distribution function that is a highly- 
irregular function of the I^s. This latter procedure might strike one as contrived and, in 
any event, one knows that, in many cases, such irregular /o's are linearly unstable. Thus, 
e.g., it is well known that, for a spherical equilibrium with fo = fo(E, J 2 ), stability or lack 
thereof often correlates with the sign of the partial derivatives dfo/dE and/or dfo/dJ 2 . In 
particular, population inversions can trigger instabilities. 

If any discrete construction based on Schwarzschild's method is to be reasonable, there 
must be a sense in which, as the discretisation of the density becomes more refined and as the 
number of building blocks becomes larger, the solution constructed numerically converges 
towards a continuous self-consistent equilibrium. However, identifying the precise sense in 
which this is so would most likely be very difficult. Mathematically, establishing such a 
convergence would probably involve a study of sequences of discrete Banach spaces, along 
the lines that have been used to study the convergence of finite difference schemes for 
solving partial differential equations. In that setting, a good deal is known about linear 
differential equations but, if one incorporates nonlinearities and/or allows for an integro- 
differential equation - recall that the collisionless Boltzmann is a quadratically nonlinear 
integro-differential equation - things become much harder! 

It is evident, both intuitively and from painful experience (cf. Siopis, Athanassoula, 
and Kandrup 1997), that it is easier to approximate comparatively smooth quantities on 
a finite lattice than quantities that manifest intricate structures on a variety of different 
scales. For this reason, one might expect that it is much easier to construct a satisfactory 
numerical representation of an fo that is a function only of smoothly varying global isolating 
integrals than an fo that depends sensitively on "local" integrals that manifest the details 
of the complex phase space structure associated with a generic nonintegrable potential. 
Moreover, even if one allows for local integrals, the numerical construction should be more 
straightforward if, for example, on a constant energy hypersurface, the phase space popu- 
lation is reasonably smooth, e.g., perhaps avoiding chaotic regions but weighting different 
regular regions in a fashion that varies smoothly with their phase space location. 

Suppose that there is in fact a "true" fo involving local integrals, generated (in principle) 
as an exact time-independent solution to the collisionless Boltzmann equation, to which one 
has constructed a latticized fo via some analogue of Schwarzschild's method, and that this 
latticized fo has been used to generate an ensemble of initial conditions to populate an N- 
body realisation of the model. There are then two potentially serious sources of error: (1) 
The latticized approximation to fo could miss important microscopic structures associated 
with local integrals. If the "true" fo is a function only of smoothly varying integrals like the 
energy E, allowing for as few as 20 different energies, as did Schwarzschild (1979), Merritt 
and Fridman (1996), and Siopis (1997), may be adequate to capture the essence of the 
analytic model. If, however, fo involves a complex combination of local integrals as well as 
the energy, allowing for 20 values may not be enough. (2) Even if most/all the important 
microscopic structures are adequately represented in the discretised model, the iV-body 
realisation could fail to sample them adequately. Even for very large particle number, 
N ~ 10 6 or more, there is no guarantee that a complex phase space will be adequately 
sampled. 

If, however, it is difficult for galactic astronomers to construct iV-body realisations of 
"complex" equilibria fo that involve local integrals in a highly irregular way, nature too 
may find it hard. If N ~ 10 6 is not enough to generate a "fair" sampling, N ~ 10 11 may 
also be inadequate to probe the complex phase space associated with a smooth potential 



12 



generated as a Boltzmann equilibrium. Both statistical fluctuations, which will obviously 
be present, and small non-Hamiltonian irregularities, which can be important in complex 
Hamiltonian systems (cf. Lichtenberg and Lieberman 1992) may tend to "smooth out" a 
complex would-be equilibrium into something substantially simpler. 
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